home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
HPAVC
/
HPAVC CD-ROM.iso
/
3DGPL.ZIP
/
3DGPL
/
TEXT
/
2.TXT
< prev
next >
Wrap
Text File
|
1995-06-22
|
46KB
|
1,117 lines
3-D transformations.
----------------------
"Several large, artificial constructions
are approaching us," ZORAC
announced after a short pause.
"The designs are not familiar, but they
are obviously the products of
intelligence. Implications: we have been
intercepted deliberately by a means
unknown, for a purpose unknown, and
transferred to a place unknown by
a form of intelligence unknown.
Apart from the unknowns, everything
is obvious."
-- James P. Hogan, "Giants Star"
3D transformations would assume the central position in the engine.
After all, this is something which would allow us to see objects from
different sides, manage the three dimensional space etc. The
simplest but by no means least important transformations are
translation and scaling. The rotation transformation would be little
bit more complicated to understand and implement but fortunately
not much. Talking of all three transformations we can either consider
them from the point of an object (collection of points in 3D)
being transformed, or a space itself being transformed. In
some cases it is more productive to think one way in some another way.
Translation and scaling.
------------------------
Translation. This transformation moves points of an object to a new
specified position pic 2.1, On the other hand we may think of it as a
system of references being moved in the opposite direction pic 2.2:
Y ^ ^
| * A'
| ^ * B' |
| | ^ Y'^
|A*--------+ | y component | |A*
| B*--------+ | B*
| x component | |
---------+-------------------> X - - - - - - | - - - - - - -> X
0| | 0|
| -----------+------------->X'
| 0'| |
| |
| |
| |
pic 2.1 pic 2.2
This way if it is a collection of points we are moving, Distances
between them would not change. Any move across 2D space can separate
into 2 components - move along X and move along Y. It would take
3 components in 3D space - along X, Y and Z. Important place we
may need translation is to , for example , move (0,0) of the screen
from upper left corner of the screen to it's physical centre. So
the translation function might look like:
void T_translation(register int *from,register int *into,int number,
int addx,int addy,int addz
)
{
register int i;
for(i=0;i<number;i++)
{
(*into++)=(*from++)+addx;
(*into++)=(*from++)+addy;
(*into++)=(*from++)+addz; /* translation */
}
}
Another transformation is scaling: blowing up of distances between
points. Or we can think of it as a contraction or blowing of the space
itself pic2.3:
A'* Y ^
\ |
\ |
* B'* *C'
A ^ /
B* *C
|
--------------+-----------------> X
0|
|
|
pic 2.3
The obvious usage of this transformation: trying to fit 2000x2000 object
into 200x200 window - evidently every point of this object would have to
have coordinates scaled to fit in the window range. In this case scaling
is just multiplication by 1/10. The object is contracted that way. On the
other hand we can think that it was actually the window space expanding to
enclose the object (we would have to abstract far enough from computer
hardware today on the market to do that, however). In general case
we can introduce separate factors for every axis, if this factor is
greater then 0 and less then 1 we would contract the object if it would
be greater then 1 the object would be expanded. The code is once again
trivial:
void T_scaling(register int *from,register int *to,int number,
float sx,float sy,float sz
)
{
register int i;
for(i=0;i<number;i++)
{
(*to++)=(int)(*from++)*sx;
(*to++)=(int)(*from++)*sy
(*to++)=(int)(*from++)*sz;
} /* scaling */
}
Planar case rotation.
---------------------
Of those three transformations rotation is by far the most complicated
one, let's first consider planar - 2-D case of rotation. there is some
simple trigonometry involved, sin and cos functions. just for those
who don't deal with this subject every day sin and cos for a triangle
on pic 2.4 are defined as:
+
/|
l/ |
/ |y
/alp|
+----+
x
pic 2.4
x y
sin(alp) = --- cos(alp) = ---
l l
And so if we do know sin(alp) or cos(alp) we can calculate segments x and y
knowing segment l:
x = l*sin(alp) y = l*cos(alp)
This can be considered as seeking projections of a point for which we
know the distance from the beginning of coordinates onto orthogonal
axes pic 2.5:
^
Y |
+----*
| /|
y| /l|
| / |
|/alp|
------+----+---->X
| x
|
pic 2.5
Now, let's consider the situation when we rotate the system of references
counterclockwise by the angle alp. What would be coordinates of point
A in the turned system of references Y'X' if it's coordinates in the
original system of references were (x,y)?
^ Y X'
| /
Y' Y/*\--*A /
\ / | \| /
\ / | |\ /
Yy'\ | | /Yx'
\ | |
\| /\ Xx'
-------------+---+-------------------> X
/| \/ X
/ Xy' \
/ | \
/ | \
/ | \
|
pic 2.6
Using the sin/cos formulas from above we can find projections of
x and y (that is of the projections of the point onto original axis's)
onto new axis's X'and Y'. Let's sum projection of x onto X' and projection
of y onto same axis's and do the same thing with projections of x and y
onto Y'. Those sums would be just what we are looking for, coordinates
of the point A in the new system of references X'Y'
i.e.
X' = Yx'+Xx' = y*sin(alp) + x*cos(alp)
Y' = Yy'+Xy' = y*cos(alp) + (-x*sin(alp))
Note the sign of Xy', just from looking on the pic2.6 one can see that
in this case x is being projected onto the negative side of Y'. That's
why the negative sign. So those are the transformation formulas for
counterclockwise rotating system of reference by an angle alp:
x'= y*sin(alp) + x*cos(alp)
y'= y*cos(alp) - x*sin(alp)
On the other hand we can consider it as point A itself rotating
around zero clockwise.
3-D rotation.
-------------
This subject is complicated only because it is a bit hard to visualize
what is happening. The idea on the other hand is simple: applying three
2-D rotations one after another so that each next works with the results
obtained on the previous stage. This way every time we are applying new
transformation we are working not with the original object but the
object already transformed by all previous rotations.
Each of 2D rotations this way is bound to some axis around which two
other axes would rotate.
There are few important considerations influencing the way we would
find the formulas: 1) What kind of system of references we have. What
is the positive direction of the rotations we would be applying; and
2) In what order we want to apply 2D rotations.
System of references:
---------------------
Well, we do have freedom of choosing directions of axes, for one
thing, we are working in orthogonal system of references where each
axes is orthogonal to the plane composed of remaining two. As
to where the positive direction of every axis point we have freedom
to decide. It is customary for example to have positive side of Y
directed upside. We don't actually have to follow customs if we really
don't want to, and remembering that the colourmap, we would be doing all
drawings into, have Y directed down we might want to choose
corresponding direction here. However having it pointing up is more
natural for us humans, so we might save time on debugging a while
afterwards. Let's point it up. As to the Z let's direct it away from
us and X to the right.
Y^ Z
| /
| / alp
/|<---+
| |/ |
-------|-+-----------> X
bet | | __
V/| /| gam
/----+
/ |
/ |
|
pic 2.7
As to the rotations, let's call the angle to turn XY plane around Z axis
as gam, ZY around X as bet and ZX around Y as gam, pic 2.7
Transformation order.
---------------------
The problem with transformation order is that the point consequently
turned by gam-bet-alp won't be necessarily at the same place with
where it might go if turned alp-bet-gam. The reason is our original
assumption: each next transformation works on already transformed
point by previous 2D rotations. That is, we are rotating coordinates
around moving axes. On the picture 2.7 we can see that rotation
bet is performed around axis x but the next rotation alp is performed
around axis z which after application of previous transformation
would move.
+ +----+
z / /| / / alp
/ + V bet \
------------+-------|---- x ---------+------------ x
/ | \
+----+ z
| / |
| V alp
pic 2.8
The implication is: we have to know with respect to what each
next rotation is performed, that is what is the order of applications
of 2D rotations. Let's think how we are coordinating the surrounding
world? First we think of the direction. Then we can tilt our head up
and down, and finally from there we can shake it from left to right.
Now, when we are tilting head we have already chosen the direction.
If we first would tilt head then the directional rotation to be
performed after would not be parallel to the ground but rather
perpendicular to the imaginary line being continuation of our
tilted neck. ( No responsibility hereby assumed for all direct
or consequential damage or injury resulted from doing that ) So
it all comes to directional rotation first - gam in our case.
pitch second - bet; roll last - alp. Let's do just that using
obtained formula for 2D rotation:
^Z ^Y' ^Y"
| | |
|<---+ |<---+ |<---+
| gam| | bet| | alp|
| | | | | |
-+----+--->X -+----+--->Z' -+----+--->X"
| | |
x'=z*sin(gam)+x*cos(gam)
y'=y
z'=z*cos(gam)-x*sin(gam)
x"=x
y"=y*cos(bet)-z*sin(bet)
z"=y*sin(bet)+z*cos(bet)
x"'=y*sin(alp)+x*cos(alp)
y"'=y*cos(alp)-x*sin(alp)
z"'=z
pic 2.9
That's basically it, using those formulas we can apply 3D rotation
to coordinates x,y,z and get x"',y"',z"'. We can see here 12
multiplication's, The question is can we reduce this number? Let's
try getting rid of temporary variables x',y',z',x",y",z". We can do it
by just substituting each occurrence of them by their real meaning:
First let's try obtaining x",y" and z" expressed via x,y,z:
y"=y'*cos(bet) - z'*sin(bet)
y"=y*cos(bet) - (z*cos(gam)-x*sin(gam)) *sin(bet)
y"=y*cos(bet) - z*cos(gam)*sin(bet) + x*sin(gam)*sin(bet)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - -
x"=x'
x"= z*sin(gam) + x*cos(gam)
- - - - - - - - - - - - - - -
z"=y'*sin(bet) + z'*cos(bet)
z"=y*sin(bet) + (z*cos(gam)-x*sin(gam)) *cos(bet)
z"=y*sin(bet) + z*cos(gam)*cos(bet) - x*sin(gam)*cos(bet)
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
Now using that we can express x"',y"' and z"' via x,y,z:
y"'=y"*cos(alp) - x"*sin(alp)
y"'=(y*cos(bet) - z*cos(gam)*sin(bet) + x*sin(gam)*sin(bet))*cos(alp) -
(z*sin(gam) + x*cos(gam)) *sin(alp))
y"'=y*cos(bet)*cos(alp)-z*cos(gam)*sin(bet)*cos(alp)+x*sin(gam)*sin(bat)*cos(alp)
-z*sin(gam)*sin(alp) - x*cos(gam)*sin(alp)
y"'=x* [ sin(gam)*sin(bet)*cos(alp) - cos(gam)*sin(alp)] +
y* [ cos(bet)*cos(alp) ] +
z* [-cos(gam)*sin(bet)*cos(alp) - sin(gam)*sin(alp)]
-----------------------------------------------------------
x"'= y"*sin(alp) + x"*cos(alp)
x"'=(y*cos(bet) - z*cos(gam)*sin(bet) + x*sin(gam)*sin(bet))*sin(alp) +
(z*sin(gam) + x*cos(gam)) *cos(alp)
x"'=y*cos(bet)*sin(alp)-z*cos(gam)*sin(bet)*sin(alp)+x*sin(gam)*sin(bet)*sin(alp)
z*sin(gam)*cos(alp) + x*cos(gam)*cos(alp)
x"'=x* [sin(gam)*sin(bet)*sin(alp) + cos(gam)*cos(alp)] +
y* [cos(bet)*sin(alp) ] +
z* [sin(gam)*cos(alp) - cos(gam)*sin(bet)*sin(alp) ]
----------------------------------------------------------
z"'=z"
z"'=x* [- sin(gam)*cos(bet)] +
y* [sin(bet)] +
z* [cos(gam)*cos(bet)]
----------------------------------------------------------------
This really appear to be more complicated that what we had before.
It appear to have much more multiplications than we had. But if we
look on those resulting formulas we can see that all coefficients
in square brackets can be calculated just once so that transformation
of a point would look like:
x"'=x*mx1 + y*my1 + z*mz1
y"'=x*mx2 + y*my2 + z*mz2
z"'=x*mx3 + y*my3 + z*mz3
It can be seen that it takes just 9 multiplication. Of course calculating
all the coefficients takes 16 multiplications:
mx1= sin(gam)*sin(bet)*sin(alp) + cos(gam)*cos(alp)
my1= cos(bet)*sin(alp)
mz1= sin(gam)*cos(alp) - cos(gam)*sin(bet)*sin(alp)
mx2= sin(gam)*sin(bet)*cos(alp) - cos(gam)*sin(alp)
my2= cos(bet)*cos(alp)
mz2= -cos(gam)*sin(bet)*cos(alp) - sin(gam)*sin(alp)
mx3= -sin(gam)*cos(bet)
my3= sin(bet)
mz3= cos(gam)*cos(bet)
But if we have 100 points to turn. Original method would take
12*100=1200 multiplications. This one would take 9*100+16=916 since
coefficients are calculated just once for all points to transform. One
more consideration having to do with optimization, in most cases it makes
sense to measure angles as integers. We don't usually care of fractions
smaller then 1 degree. So we don't really need sin and cos for all the
real number range. That's why we can calculate sin/cos just once for
all 360 degrees before any rotations are performed and put the obtained
values into arrays so that next time we would need any of them we
wouldn't have to call a functions which probably uses a power series
with few multiplications to calculate each. (Processors nowadays come
with floating point units which can calculate sin/cos pretty fast
but unlikely faster then single array access (this might become false
though). By the way we don't really need 360 degrees. This number is
not very convenient. We can go with full angle divided into 256 pseudo
degrees. The reason we would need just one unsigned char to store the
angle and whenever we go beyond 255 we simple wrap around to zero,
saving this way a conditional statement we would need otherwise in
the case of 360 degrees.
#define T_RADS 40.7436642 /* pseudo grads into rads */
float T_mx1,T_mx2,T_mx3,T_my1,T_my2,T_my3,T_mz1,T_mz2,T_mz3;
float T_sin[256],T_cos[256]; /* precalculated */
void T_init_math(void)
{
int i;
for(i=0;i<256;i++)
{
T_sin[i]=sin(i/T_RADS);
T_cos[i]=cos(i/T_RADS);
}
}
Now finally to the code to perform 3D rotation. First function is building
set of rotation coefficients. Note that T_mx1, T_mx2 etc. are global
variables. The other function: T_rotation uses those coefficients, so
they better be ready when later is called.
void T_set_rotation_gam_bet_alp(int alp,int bet,int gam)
{
T_cosalp=T_cos[alp];
T_sinalp=T_sin[alp];
T_cosbet=T_cos[bet];
T_sinbet=T_sin[bet];
T_cosgam=T_cos[gam];
T_singam=T_sin[gam]; /* initializing */
T_mx1=T_singam*T_sinbet*T_sinalp + T_cosgam*T_cosalp;
T_my1=T_cosbet*T_sinalp;
T_mz1=T_singam*T_cosalp - T_cosgam*T_sinbet*T_sinalp;
T_mx2=T_singam*T_sinbet*T_cosalp - T_cosgam*T_sinalp;
T_my2=T_cosbet*T_cosalp;
T_mz2=-T_cosgam*T_sinbet*T_cosalp - T_singam*T_sinalp;
T_mx3=-T_singam*T_cosbet;
T_my3=T_sinbet;
T_mz3=T_cosgam*T_cosbet; /* calculating the matrix */
}
void T_rotation(int *from,register int *to,int number)
{
register int i;
register int xt,yt,zt;
for(i=0;i<number;i++)
{
xt=*from++;
yt=*from++;
zt=*from++;
*to++=(int)(T_mx1*xt+T_my1*yt+T_mz1*zt);
*to++=(int)(T_mx2*xt+T_my2*yt+T_mz2*zt);
*to++=(int)(T_mx3*xt+T_my3*yt+T_mz3*zt); /* performing the operation */
}
}
Representing transformations in matrix form.
--------------------------------------------
Matrix is a table of values:
| a b c |
A=| d e f |
| g i j |
Special cases are column and row vectors:
| x |
V=| a b c | C=| y |
| z |
Row can by multiplied by a column producing a value, a scalar:
| x |
| a b c |*| y | = a*x + b*y + c*z
| z |
Matrix can be multiplied by another matrix: each value of the
result was obtained by multiplying respective rows and columns
of the arguments (since columns and vectors are to have the same
dimension in order for multiplication on them to be defined it
automatically places restrictions on dimensions of matrices that
can be multiplied):
| x | | k |
| | a b c |*| y | | a b c |*| l | |
| | z | | m | |
| |
| a b c | | x k | | | x | | k | |
| d e f |*| y l | = | | d e f |*| y | | d e f |*| l | |
| g i j | | z m | | | z | | m | |
| |
| | x | | k | |
| | g i j |*| y | | g i j |*| l | |
| z | | m |
We may be interested in multiplying a row vector by a 3x3 matrix:
| a b c |
| x y z |*| d e f | = | x*a+y*d+z*g x*b+y*e+z*i x*c+y*f+z*j |
| g i j |
Looks familiar? Yes, the rotation coefficients we've calculated
lie nicely in the 3x3 matrix and so we can think of 3D rotation
as of matrix multiplication of "rotation" matrix onto original
vector producing result vector. We can also try fitting translation
and scaling into a matrix, to do that for translation, however, we
would have to increase dimensions of the coefficient matrices and
the argument vector:
Translation:
| 1 0 0 0 |
| x y z 1 |*| 0 1 0 0 | = | x+tx y+ty z+tz 1 |
| 0 0 1 0 |
| tx ty tz 1 |
Scaling:
| sx 0 0 0 |
| x y z 1 |*| 0 sy 0 0 | = | x*sx y*sy z*sz 1 |
| 0 0 sz 0 |
| 0 0 0 1 |
The advantage of the matrix approach is that if we have several
transformations each represented in matrix form:
| X'| = (| X |*| A |) *| B |
Where A, and B are transformation matrixes and | X | the argument
vector to be transformed we can:
| X'| = | X |* ( | A |*| B | )
calculate the "concatenated" matrix | K | = | A |*| B | and each
following coordinate just transform as:
| X'| = | X |*| K |
By the way by representing each of 2-D rotations in 3x3 matrix form we
can get the formulas we have found for 3-D rotation by calculating a
concatenated matrix.
So, matrix approaches might be very effective in instances where
there are a lot of consecutive transformation. For example if we
have consecutive rotations it makes sense to calculate their
concatenated matrix and use it from then on. It would save
us quite a few multiplication, one verses several consecutive matrix
multiplication for each point to transform.
And in the run-time it is preferably to use pre-calculated
formulas for concatenated matrix coefficients then trying to get it
every time by multiplying respective transformational matrices.
The reason 3 consecutive matrix multiplication of, for instance, 3x3
matrices take 9+9+9=27 multiplication's. But since lot of values in
those matrices are 0, it is taken into account in pre-calculated
coefficient formulas which allow to produce the very same results
with just 16 multiplications.
As to the translation transformation matrix and especially increasing
the dimension from 3 to 4 it looks much less natural. 4x4 matrix means
more multiplications too, of course some of the value in 4x4 would
always be zero so we can optimize (cripple) the routine to original 9
multiplications and 3 additions. In one word, there's a freedom
to decide what approach to take.
perspective.
------------
I believe, it is not hard to visualize what perspective is, perhaps the
first association lot of us would have, is this empty straight street
with identical buildings on both sides disappearing in infinity. ( At
least people residing in newer neighborhoods of Moscow or Toronto know
what I am talking about ) Why do we need this transformation? Just to
achieve landscape realism, after all this same street would look pretty
weird if the road would not collapse into a dot and buildings further
away from us won't look smaller. So weird it might be even hard to
picture.
Perspective might not be needed for CAD style applications, but even
there it would, most likely, be implemented in one viewing mode or
another.
Before we can implement this transformation it is a good idea to understand
physics of what is happening. A ray from some point in space is being
caught by the viewer's eye. This ray had intersected a plane located
in front of the viewer. If we can find an intersection point and plot the
dot there, for the viewer it would be the same as if the ray from the
plotted dot was actually coming from the initial point in space.
* A
| /
A'| /
*
/ | B'
*---*-----------------------* B
\ |
*\
C'| \
| \
* C
pic 2.10
Rays from different points in space , in this model, all converge
in one point - viewer's eye. The math here is quite straightforward
pic 2.11:
screen plane->| *
| /| A (X,Z)
| / |
|/ |
^ * |
| /|X' |
/ | |
|/ | |
*---+---+- - - - ->
0 | Z
| |
focus distance ->|---|<-
pic 2.11
Let's consider 2-D case first: the viewer's eye is located at point 0
the distance between viewer's eye and the screen would be called -
"focus". The problem is to determine at what point the ray going from
point A into the viewer's eye would intersect the screen. We would have
to plot the dot just there. This is very simple, just yet another
case of solving similar triangles. Just from the fact that the angle
at 0 is the same for both triangles, and if two angles are the same their
tangents would be too. We have:
X' X X*focus
----- = ----- => X'= -----------
focus Z Z
Same considerations apply for Y dimension, together describing 3-D case:
Y'=Y*focus/Z X'=X*focus/Z
It is a bit of a mystery what is happening with Z coordinate. Basically
after the perspective transformation we would be performing rendering
onto the screen which is 2-D surface we would need X and Y but hardly
Z, so we might as well have Z'=0. However when trying to render multiface
objects we would need to know the depth (Z) of all polygons so that
we can decide which are visible and which are obscured. So we might
have Z'=Z that is just leave old value. Then again by doing transformation
on X and Y and leaving Z unchanged we actually may allow for depth
relation to change. That is a polygon obscuring another one in the original
space, before projection, can actually come out behind or partly behind
with this kind of transformation pic 2.21.
^ Z
same | * *
deapth | * / * /
before | A / / //
and | * / /* A
after | /B / B
| * *
| before after
pic 2.12
To avoid that we can apply some non linear transformation on Z,
like Z'=C1+C2/Z, where C1, C2 are some constants.
In the enclosed code only X and Y are being transformed. The implementation
of all those things is that easy I won't even put a code example here.
However, there are few possible questions: first: what should be the value
for "focus distance" ?
\ / \ /
\ / \ /
\ / \ /
---I\--------/I--- ---I\--------/I---
\ / \ /
\/ \ /
\ /
\/
pic 2.13
pic 2.13 shows it's physical sense, focus basically determines the view
angle. when it is smaller - angle is wider, when it is bigger angle
is smaller. It is helpful to think of this distance as being measured
in screen pixels, this way for 320 pixel display if focus is chosen
to be 180 the view angle we are getting is 90'. Perspective
transformation might produce at first, weird looking results, with
distortions sometimes similar to that of wide-range photographic
pictures. One has to experiment a bit to find what's looking best,
values giving view angle somewhere between 75-85' would probably
look not that terrible, but that depends of scene and screen geometries.
clipping problem...
-------------------
Where transforms point with Z equal to 0? since we are dividing by
Z that would be infinity, or at least something producing divide error,
in this computer-down-to-earth-world. The only way to deal with this
is to guarantee that there would be no coordinates with such Zs. another
thing, calculations for points with negative Z would produce negative
coordinates, what we would see is objects (or even parts thereof) flipped
over, but then again objects with negative coordinates are behind the
viewing plane and effectively behind the viewer, so this is something
we can't see and better make sure this is not being rendered. The way
to do it is by applying 3-D clipping.
The perspective transformation can also be represented as a matrix,
however, let's think of all transformations discussed up to this moment
a point would have to go through in order to be rendered. First we would
rotate the world with respect to viewer orientation, then before we can
apply perspective transformation we would have to do clipping and
at least get rid of the vertices behind the viewing plane, and only then
can we apply perspective transformation. Representing clipping in
matrix form is not the most pleasant thing there is. And we've already
discussed that matrices would be of use when we have several consecutive
transformations so that all of them would be represented by one matrix.
In this case however rotation and perspective are separated by the clipping.
If we are sure that none of the coordinates would ever get behind the
viewing plane, we don't actually need clipping, that's the instance to
represent everything in the matrix form, otherwise it might not be
a good idea neither in terms of code complicity nor performance.
fixed point math.
-----------------
In previous chapters in several places we didn't have another choice but to
use floating point multiplications. (sin and cos are after all real-valued
functions) However this is still fairly expensive operation, integer
multiplications and divisions would usually be faster. However since we
have to represent numbers with fractions it appeared we didn't have another
choice but to use floating point operations. The question is, can we
represent fractions using just integers? One might think there's a
contradiction in the question itself, however we would see that, as in
many other cases, we can use one mechanism to represent something else.
This is the case here, We may use fixed point numbers representing them
as integers, and with small amendments we can use regular integer additions
multiplications etc, to perform fixed point operations. First of all how do
we represent integers?
+------+------+------+------+------+
| 1 | 0 | 0 | 0 | 1 |
+------+------+------+------+------+
4 3 2 1 0
2 =16 2 =8 2 =4 2 =2 2 =1
pic 2.14
With any counting system, digits in multi digit numbers would have
different weight, so to speak. With decimal numbers this weight
would be some power of ten, that is 102 is actually:
10**2*1 + 10**1*0 + 10**0*2 == 100+0+2 ==102
Same with binary numbers, the only difference, weight would be
some power of 2. this way number on pic 2.13 would be:
2**4*1 + 2**3*0 + 2**2*0 + 2**1*0 + 2**0*1 == 16+0+0+0+1 == 17
Now, it is quite natural for us to place a decimal point into a decimal
number, placing a binary point into a binary number should not
seem more difficult. How do we interpret the decimal point? It is
just that for the numbers to the right of the decimal point we
pick negative power of ten. That is:
102.15 == 10**2*1 + 10**1*0 + 10**0*2 + 10**-1*1 + 10**-2*5 ==
1 5 15
== 100 + 0 + 2 + ---- + ---- == 102 + -----
10 100 100
In this representation the point separates negative powers from
zero power and unlike numbers represented in exponential form - the
floating point numbers, in this case the point never changes it's
position.
Nothing should stop us from extending the same scheme onto binary
numbers:
+------+------+------+------+------+
| 1 | 0 | 0 | 1 | 1 |
+------+------+------+------+------+
2 1 0 -1 -2
2 =4 2 =2 2 =1 2 =1/2 2 =1/4
pic 2.15
Using the very same technique we can see that number represented
on the pic 2.15 can be considered also as:
2**2*1 + 2**1*0 + 2**0*0 + 2**-1*1 + 2**-2*1 == 4+0+0+1/2+1/4
== 4 + 3/4
Numbers in what range can we represent? Resorting once again to
analogy with decimal numbers it can be seen that two digits
to the right from the decimal point can cover the range of 0.99
in 0.01 steps. Numbers smaller then minimal 0.01 precision step
just can't be represented without increasing number of digits
after the decimal point. Very same thing is happening with binaries
in the example at pic 2.15 since we have two binary digits- minimal
number we can represent is 1/4 ( 0.01 (bin) == 1/4 (dec) ) and
again numbers with higher precision can be represented only by
increasing number of binary digits after the binary point.
Now to the representation of negative numbers, there actually
several ways to do that but effectively today one has prevailed,
this is so called 2's compliment form. The idea is to negate
the weight of the leftmost digit pic 2.16:
+------+------+------+------+------+
| 1 | 1 | 1 | 1 | 1 |
+------+------+------+------+------+
4 3 2 1 0
-2 =-16 2 =8 2 =4 2 =2 2 =1
pic 2.16
the number represented above is considered to be:
-2**4*1 + 2**3*1 + 2**2*1 + 2**1*1 + 2**0*1== -16+8+4+2+1 = -1
The advantage of this approach, we don't have to change addition
and subtraction algorithm working for positive integers in order
to accommodate 2's compliment numbers. Why is that? just by the
virtue of natural overwrap of integer numbers the way they are
represented in the computers. If we add 1 to the maximum number
we can represent we would obtain a carry from the most significant
(leftmost) bit, which is ignored, and the value of exactly 0.
-1 in signed representation is exactly the maximum possible value
to represent in unsigned representation, -1+1==0 indeed. (One can
check this to be true for all numbers and by more formal means, but
since this is not the point I would ignore it) Again, addition and
subtraction don't have to be changed to accommodate 2's compliment
negative numbers (multiplication however should, that's why there
are instructions for both signed and unsigned multiplications
in most computers)
Since the leftmost digit in 2's compliment representation carries
negative weight and that's exactly the one with highest power the
minimum negative number possible to represent in the example
above would be 10000 (bin) = -16 all the other digits have
positive weights so maximum possible positive number would be
01111 (bin) = 15 but this slight asymmetry is not a problem in
the majority of cases.
However, values of sin and cosin functions are between 1 and -1 so
to represent them we might choose format with just one integer
field:
+------+------+------+------+------+
| 0 | 1 | 1 | 1 | 1 |
+------+------+------+------+------+
0 -1 -2 -3 -4
-2 =-1 2 =1/2 2 =1/4 2 =1/8 2 =1/16
pic 2.17
but, because of the asymmetry described above there would be a
value for -1 (1000) but there would be no pure 1, just it's
approximation (01111): (pic 2.17) 1/2+1/4+1/8+1/16, Then again, for
most graphics applications when we have, say, 15 fields representing
fractions this would be close enough and won't cause a lot of problems.
As you can see, the very same rules work for 2's compliment whether
with fixed point or without.
multiplication of fixed point numbers.
--------------------------------------
Let's consider what is happening when we are multiplying two
decimal numbers, for example we have to multiply an integer
by a number with a decimal point and we would need just
an integer as a result:
1.52
* 11
-----
1 52
+ 15 2
-----
16.72
Actual result of multiplications has as many digits after
the decimal point as there were in both arguments. Since
we need just an integer as a result we would discard numbers
after the point effectively shift digits to the right by
2 in this case. We can do the same with fixed point binary
numbers, it means that we have to readjust precision after
the multiplication. Say, if we are multiplying an integer
by a number having 8bit considered to be after the binary
point the result would also 8bit after the point. If it
is just the integer part we are interested in, we have to
shift the result right 8 bit destroying all fractional
information. Similar things are happening with division
except that if we want to divide two integers and obtain
a fixed point result the argument to be divided has to
be shifted left- added fixed point fields, so that effectively
we are dividing a fixed point number by an integer.
How to implement all that? Well, first of all we have to
decide what kind of values and what operations on them
we would need. Sometimes it is nice just to have general
purpose fixed point math library, on the other hand we
may choose several different formats of fixed point numbers
one for every special place. Another thing we have to
consider is: would we be able to represent all the operations
using just high level C constructs or it would be necessary
to descend into assembly instructions. The reason we might
consider later is: in most assemblies the result of
multiplication has twice as many bits as each of arguments,
and since we would occupy some bits with fractions we really
need all bits we can get. Another thing, the result of
multiplication would often be located in two registers, So
we can organize operations the way so that instead of adjusting
result by shifting we would just take it from the register
carrying higher bits, -effectively doing zero-cost shifting
right by the bit length of the register. On the other hand,
do we really want to compromise code portability by coding
in assembly? Obviously, there are valid reasons to answer that
both ways. In particular case of this text and associated
code we are interested in good portability so let's try to
stay within straight C boundaries.
Where would we be using fixed point? first of all - in 3-D
transformations, to keep rotation matrix coefficients. We
would need two kinds of multiplication: In calculation of
the coefficients: fixed * fixed producing same precision
fixed and in actual coordinate transformation fixed * int
producing int. lets define fixed to be:
typedef signed_32_bit fixed; /* better be 32bit machine */
#define T_P 10 /* fixed math precision */
the T_P would be the number of bits we are allocating to
keep fractions. Let's make a little function to transform the
floating point values in the range [-1,1] into fixed point
values:
fixed TI_float2fixed(float a)
{
fixed res=0;
int i,dv,sign;
if((sign=a<0.0)==1) a=-a;
for(i=1,dv=2;i<=T_P;i++,dv*=2)
if(a>=1.0/dv) { a-=1.0/dv; res|=(0x1<<(T_P-i)); };
if(sign==1) res=-res;
return(res);
}
In the function above we are iteratively subtracting powers
of two from the argument and using the result to set the bits
in the fixed point number. We would use this function to
precalculate the array of sin/cos values.
When calculating the rotational coefficients we would multiply
one fixed by another and obtain same precision fixed result.
again, actual result would have twice as many fractional bits
as each of the argument so same precision result would be
(a*b)>>T_P where 'a' and 'b' fixed point values and T_P number
of bits after the binary point in each. This way code would
look something like:
...
...
T_wx1=((singam*((sinbet*sinalp)>>T_P))>>T_P) + ((cosgam*cosalp)>>T_P);
T_wy1=((cosbet*sinalp)>>T_P);
T_wz1=((singam*cosalp)>>T_P) - ((cosgam*((sinbet*sinalp)>>T_P))>>T_P);
...
...
In rotation function ints would have to be multiplied by a fixed
coefficient, the result would have T_P fractional bits, since we
are interested in just integer part same rule as above applies:
...
...
*to++=((T_wx1*xt)>>T_P) + ((T_wy1*yt)>>T_P) + ((T_wz1*zt)>>T_P);
*to++=((T_wx2*xt)>>T_P) + ((T_wy2*yt)>>T_P) + ((T_wz2*zt)>>T_P);
*to++=((T_wx3*xt)>>T_P) + ((T_wy3*yt)>>T_P) + ((T_wz3*zt)>>T_P);
...
...
That's about it. The only thing, we really have to be careful with
the range of numbers the operations would work for. If we are working
with 32bit numbers and would choose 16 of them for fractional part
and one for sign bit then only 15bit would carry the integer part.
after multiplication of two fixed numbers since result has as many
fractional bits as in both arguments all 32 bits would actually carry
fractions. The integer part would be gone in the dark realm of left-carry,
and there's no way in standard straight C to get it out of there.
(when using assembly we don't have this problem since result
of multiplication has twice bits of each argument, so we really can
get trough to the integer part, however, most C compilers define the
result of multiplication to be of the same bit size as arguments
int*int==int, some compilers try to be smarter then that but we really
have to consider worst-case scenario)
And finally is it really true that integer multiplication plus shifts
is faster then floating multiplication? What follows are results of
some very approximate tests:
sparc: floating faster fixed about 1.3 times
68040: floating faster fixed about 1.5 times
rs4000: floating faster fixed about 1.1 times
rs6000: fixed faster floating about 1.1 times
i386sx: fixed faster floating about 5.1 times
floating in the test was: float a,b; a*b;
fixed in the test was: int a,b; (a*b)>>15;
As one can see processors with fast floating point units nowadays
do floating operations faster then integer once. Then again in terms
of portability integer multiplication would always be fast
enough whereas floating multiplication especially on processors
without FPUs (i386sx) might get really slow. The decision whether
to implement fixed point, this way, might get a bit tough to make.
In the provided library's source transformations are implemented
both ways but with the same function prototypes, so that one can
decide which implementation is more suitable to link in.
* * *
